Narrow depression in the density of states at the Dirac point in disordered graphene 
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The electronic properties of non-interacting particles moving on a two-dimensional bricklayer 
lattice are investigated numerically. In particular, the influence of disorder in form of a spatially 
varying random magnetic flux is studied. In addition, a strong perpendicular constant magnetic 
field B is considered. The density of states p(E) goes to zero for E — > as in the ordered system, 
but with a much steeper slope. This happens for both cases: at the Dirac point for B = and at 
the center of the central Landau band for finite B. Close to the Dirac point, the dependence of p(E) 
on the system size, on the disorder strength, and on the constant magnetic flux density is analyzed 
and fitted to an analytical expression proposed previously in connection with the thermal quantum 
Hall effect. Additional short-range on-site disorder completely replenishes the indentation in the 
density of states at the Dirac point. 

PACS numbers: 73.23.-b, 71.30.+h, 73.22.-f 



I. INTRODUCTION 

Despite enormous efforts in recent years, the elec- 
tronic properties of graphene near the neutrality point 
that separates conduction and valence bands are still un- 
der intense investigations. Theories aimed at describing 
graphene via a Dirac-likc equation have elucidated many 
intriguing effects like the physics of massless relativis- 
tic particles and Klein tunneling ^ 2 ) 3 ) 4 ' 5 !, which has been 
experimentally observed recentlyi&I The Dirac-fermion 
approach, which represents an approximation to the true 
lattice situation, is believed to be valid especially close 
to energy zero, the so called Dirac points, where the con- 
duction and valence bands touch each other and the dis- 
persion is linear. Furthermore, nearest neighbor tight- 
binding descriptions ) ' ' which particularly emphasize 
the hexagonal lattice structure of the carbon sheet, have 
proven to be extremely helpful in understanding the basic 
transport properties of this promising new material. 

Based on the observation of a peculiar quantum Hall 
effec t 11 ' 12 : 13 , it is generally accepted by now that some 
kind of disorder must be present in the experimental 
setup. The disorder influences the charge transport 
through the graphene sheet and affects the measurable 
quantities at least quantitatively^ Yet, which type of 
disorder is encountered in real samples is still completely 
unclear or only partly known in some special cases. This 
lack of knowledge is particularly unfortunate as the def- 
inite type of disorder entirely determines the physical 



properties^ 
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e.g., leading to complete An- 



derson localization in the case of short range electrostatic 
scattering potentials via chiral symmetry breaking and 
scattering between valleys i 22 : 23 For disordered systems, 
even the single particle density of states (DOS) near the 
Dirac point remains still under debate. Depending on 
the disorder type and approach, a vanishing, a finite, or 
an infinite DOS at the Dirac point has been suggested 
for graphene or related models 

Also, the interpretation of experimental results is ham- 
pered by the uncertainty regarding the precise form of 
the DOS. Recently, following earlier experimental inves- 



tigations of the Landau level splitting in high magnetic 
fields ; 35 : 36 the opening of a spin (Zeeman) gap in the 
density of states at the Dirac point has been suggested 
in the interpretation of magneto-transport measurements 
on graphene sheets^ 7 - Only if a gap, separating electron 
and hole states at the Dirac point, was assumed, the 
experimental data could be accounted for. Another un- 
explained experimental observation to be found near the 
Dirac point in the presence of a strong magnetic field is 
the divergent resistance ; 38 : 39 which has been attracting 
considerable attention lately. 

Theoretically, the opening of a mobility gap within 
the central Landau band has been recently discov- 
ered by means of detailed two-terminal conductance 
calculations^ It was found that with increasing disorder, 
the critical energies where the plateau transitions of the 
Hall conductivity take place, move apart. This splitting 
of the central conductance peak unveiled the existence 
of an extra chiral quantum phase transition occurring at 
zero energy with critical properties that differ from those 
of the quantum Hall transitions^ 

In the present work, the single particle density of states 
is calculated numerically for the same bricklayer lattice 
model. The presumed ripple disorder is modeled by a 
spatially varying random magnetic flux with zero mean, 
pointing perpendicular to the two-dimensional lattice. 
Also, an additional constant magnetic field is applied 
that leads to the formation of Landau bands. It is shown 
that in the disordered case, the density of states goes to 
zero at the Dirac point not only in the absence of a per- 
pendicular magnetic field. Rather, a narrow suppression 
in the DOS is obtained also in the presence of a finite 
magnetic field within the lowest (central) Landau band. 
This unexpected feature depends essentially on the dis- 
order strength, on the system size, and on the strength 
of the perpendicular constant magnetic field. Due to the 
neglect of electron spin in the model Hamiltonian, this 
outcome can not be attributed to a Zeeman splitting, 
but must originate from chirality and a disorder-induced 
interaction between the two sub-lattices. This consider- 
ation is confirmed by the observation that the addition 
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FIG. 1: (Color online) The ensemble averaged density of 
states of a two-dimensional bricklayer lattice with random 
flux disorder strength / = 0.05, 0.2, 0.5, and 1.0. With in- 
creasing /, the van Hove singularities disappear. The size of 
the bricklayer system is L x x L y — 64 x 128 a 2 . 



of short-range potential disorder completely destroys the 
DOS-depression near the Dirac point. 



II. BRICKLAYER MODEL 

In the present study, the two-dimensional honeycomb 
lattice responsible for the peculiar electronic properties 
of graphene is replaced by a bricklayer model^^ which 
shares the same topology as the hexagonal lattice. The 
bricklayer lattice is bi-partite and consists of two sub- 
lattices that can be constructed by rectangular unit cells 
of size 2a x a placed along the rr-dircction. The unit 
cell contains two sites connected by a bond of length a. 
Each site on one sub-lattice is attached to three neigh- 
bors belonging to the other sub-lattice by two bonds in 
the ±a;-direction and one alternating bond in the in- 
direction. 

A tight-binding Hamiltonian for non-interacting parti- 
cles with nearest neighbor transfer energy V in the pres- 
ence of perpendicular magnetic fields is defined by 
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FIG. 2: (Color online) The energy dependence of the den- 
sity of states near the Dirac point with random flux disorder 
strength f/(h/e) = 0.2, 0.3, 0.4, 0.5, 0.6, 0.75, and 1.0. The 
steepest DOS tails belong to f /(h/e) = 1.0. The system size 
is L x /a = Ly/a = 96. 



arc defined by the total magnetic flux Q x y = p/q (h/e) + 
(j) x ,y threading a given plaquctte with upper left corner 
at site (x,y). The constant magnetic field is given by 
the fraction p/q of a flux quantum h/e with mutually 
prime integers p and q so that the tight-binding band 
splits exactly into 2q sub-bands, and tf> x y is the random 
flux part. The latter incorporates the effect of inhomo- 
geneous magnetic fields and mimics the disorder due to 
corrugations and ripple a 40 i 41 ' 42 present in real graphene 
sheets. In contrast to diagonal disorder, it preserves the 
chiral symmetry and ensures a finite conductivity at the 
Dirac point. 

The random fluxes are drawn from a box distribution 
~~ //2 < <j)x,y < //2 with zero mean and disorder strength 
< / /(h/e) < 1. Periodic boundary conditions are ap- 
plied in both directions to avoid edge and corner effects 
and the system size was chosen to be commensurate with 
the spatially constant magnetic field. The eigenvalues 
Ei (n) of the Hamiltonian ([I]) were obtained by direct di- 
agonalization of the N r disorder realizations and used for 
the calculation of the ensemble averaged density of states 
p(E) within an energy interval AE 

E+AE 
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and c x , y are creation and annihilation opera- 
tors of a particle at site (x, y), respectively. The prime at 
the first sum in JT]) indicates that only transfers along the 
non-zero vertical bonds are included. The second sum 
describes the movement in the horizontal chains. The 
phases, which are chosen to be only associated with the 
vertical bonds in the y-direction, 
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III. DENSITY OF STATES 

A. Magnetic field B = 

Starting with the case where the constant part of the 
magnetic flux density is zero and only the random dis- 
order part is present, the ensemble averaged density of 
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FIG. 3: (Color online) The disorder dependence of the fitting 
parameter g s (f,L) with L/a = 96 for f/(h/e) = 0.1 and 
the seven disorder values / shown in Fig. [2] On this double- 
log plot the straight line is given by g B (f) = 1580 f 2 ' 5 . In 
addition, data for L/a = 192 are shown with a power-law fit 
ffa (/)= 8702/ 2 ' 5 (.). 



states for bricklayer systems of size L x x L y = 64 x 128 a 2 
is shown in Fig. [T] for different disorder strengths /. 
With increasing /, the sharp van Hove singularities at 
E/V = ±1.0 get rounded and finally disappear. In the 
same way, the fluctuations, which can be seen for the 
smallest disorder f /{h/e) = 0.05, vanish. The latter are 
due to finite size effects. Based on an energy resolution 
of 0.02V as used in Fig. []} the other curves do neither 
depend on the shape of the system nor on the size which 
has been checked within the range 32 < L/a < 192. The 
main consequence of the increasing disorder is seemingly 
the filling of the valley in the density of states with a 
strong increase at the Dirac point E/V = 0. However, 
a closer inspection of the energy range near the Dirac 
point reveals a completely different behavior. As shown 
in Fig. O independent of disorder strength, the DOS al- 
ways goes down to zero at E/V — 0. For small random 
flux disorder, the DOS vanishes with a slope that finally 
becomes 2/(9irV 2 a 2 ) in the clean limit. With increas- 
ing disorder strength /, this slope becomes steeper and 
steeper. Since f /(h/e) = 1.0 is the strongest random flux 
disorder possible, there will always be an energy region 
around the Dirac point where the density of states goes 
to zero at E/V = 0. 

These results have been obtained with different diago- 
nalization methods including a Lanczos algorithm as well 
as standard LAPACK routines. The number of realiza- 
tions exceeded 10 4 for each disorder / and the DOS bin- 
width AE/V was 4 x 10~ 6 . 

A disorder dependent vanishing of the density of 
states was previously reported for massless random Dirac 
fermions on a two-dimensional square lattice 4 ^ and on 
a honeycomb lattice^ where non-diagonal disorder was 
introduced by real random hopping terms. These model 
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FIG. 4: (Color online) The energy dependence of the density 
of states for disorder / = h/e and system size L/a =64, 96, 
128, and 192. With increasing size, the energy range of the 
DOS-depression becomes narrower. 

systems preserve time reversal symmetry and therefore 
belong to the (chiral) orthogonal universality class. Also, 
the energy range where the DOS drops to zero is consid- 
erably broader compared with what has been found in 
our random flux bricklayer model, which belongs to the 
chiral unitary symmetry class. Furthermore, the singular 
peak observed at E = in Rcf . [34| for both the disordered 
random fermions with either random hopping or random 
gauge fields is absent in the present bricklayer situation. 

As seen from the thin black lines used to fit the curves 
in Fig. [21 neither a simple linear energy dependence, as 
observed in the clean system, nor a power-law form used 
in Refs. H3, [H andliais adequate. An additional loga- 
rithmic term similar as in the case of the class-D thermal 
quantum Hall effcc t 44 i 45 or for Dirac fermions on a hon- 
eycomb lattice with weak diagonal and bond disorder- 
seems to be an appropriate empirical function. Although 
bond and random flux disorder are different, both disor- 
der types maintain the chiral symmetry of the system. 
Therefore, we try to use the ansatz 

5^ ( 1 + ^ i > ta IS7F[)' < 4 » 

with a disorder and size dependent fitting function 
g s (f,L). The latter grows with both increasing disor- 
der strength and system size L. For square samples of 
fixed size L — (96 a) 2 as used in Fig. [5J a power-law de- 
pendence on the disorder strength g s (f, L/a = 96) oc / 
is found in the range 0.1 < f /(h/e) < 1.0, as is shown 
in Fig. [3] The same behavior is obtained for larger sizes 
L/a =128 and 192. The latter data are also included in 
Fig. El 

The overall shape of the density of states, as plotted in 
Fig. [TJ seems to be independent of the system size with 
the exception of noticeable small finite size fluctuations 
occurring only for the smallest disorder. In striking con- 
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FIG. 5: (Color online) The size dependence of the fitting 
parameter g s (f,L) for square samples and two disorder val- 
ues / = 1.0 h/e (+) and / = 0.5 h/e (o). On a double-log 
scale, the straight lines are given by g s (f /(h/e) = 1.0, L) — 
0.123 (L/a) 2 and g s (f/(h/e) = 0.5, L) = 0.0176 (L/a) 2 - 15 , re- 
spectively. 

trast, the depression of the DOS near E/V = shows a 
strong length dependence. In Fig. |4] the averaged density 
of states for disorder f = h/e is plotted for sample sizes 
L/a =64, 96, 128, and 192. Applying the function Q, 
the size dependence of g s close to E/V = can be ob- 
tained. This is shown in Fig. [5] where the size dependence 
of the fitting parameter g s (f,L) with 64 < L/a < 512 
is shown on a double-log plot for two disorder strengths 
/ = 1.0 and / = 0.5, respectively. In both cases, a power- 
law relation g s (f,L) cx L K is obtained with an exponent 
K = 2.0 for '/ = 1.0 h/e and k = 2.15 if / = 0.5 h/e. 
Because of the uncertainties due to the limited number 
of realizations, particularly for larger system sizes, and 
the restricted range 64 < L/a < 512, it is not possible to 
rule out that both exponents are the same in the limit 
L —> oo. The energy range where the function can be 
fitted to the numerical curves decreases with increasing 
size I in a similar manner as with increasing disorder 
strength /, which can be seen in FigO 

With the disorder and size dependence as identified 
above from the numerical data, our ansatz for the den- 
sity of states ((3]) used for all curves shown in Fig. [2J is 
the fitting function g s (f,L) = g c [f /(h/e)} 5 / 2 (L/a) 2 with 
only one adjustable constant g c = 0.0925. 



B. Finite magnetic field 

For a continuum Dirac model in the presence of a fi- 
nite magnetic field, the energy spectrum of the charge 
carriers is arranged into degenerate Landau levels. The 
energetically lowest Landau level appears at the charge 
neutrality point at E = 0.— Hence, instead of the den- 
sity of states going to zero, a finite DOS arises at the 
Dirac point in the disorder free system for B ^ 0. To 
account for an external perpendicular B-ficld in our lat- 
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FIG. 6: (Color online) The disorder broadening of the lowest 
(central) Landau band for square samples of size L x = L y = 
128a and constant magnetic flux density p/q — l/32(h/e)a~ 2 . 
The disorder strength is f/(h/e) = 0.02,0.03,0.05, and 0.07, 
respectively. 



tice model, a spatially constant magnetic flux is applied 
in addition to the random magnetic flux disorder. This 
type of disorder causes the central Landau level at the 
Dirac point to broaden only a littl o 18 ' 21 compared with 
the broadening of the higher Landau bands. However, 
one has to keep in mind that due to the lattice struc- 
ture, the sub-bands exhibit Harper's broadening already 
in the disorder-free system. This intrinsic broadening 
is small and disappears with decreasing magnetic field. 
The disorder broadening of the narrow central Landau 
band, which is proportional to f \fpjq~, is seen in Fig. [5] 
for p/q = 1/32 and L x = L y = 128 a. There is also 
an additional narrow structure discernible around E = 
and this was previously attributed to originate from the 
chiral critical eigenstatesJ^ 

In order to scrutinize this special feature, more than 
10 4 disorder realizations were calculated. Thereby, a bin- 
width of about 5 x 10~ 7 becomes possible leading to an 
enhanced energy resolution. As a result of this effort, one 
can see that the density of states in the center of the low- 
est Landau band is by no means constant, but is in fact 
dominated by a narrow depression (see Fig. [7]) , which de- 
pends on disorder strength, on the constant part of the 
magnetic flux density, and on the system size. This is 
an unexpected outcome and was not identified in previ- 
ous work. It would be very interesting to see whether 
or not a similar depression with a density of states going 
to zero at the Dirac point develops also in a continuum 
Dirac equation approach. The latter method is generally 
employed in graphene studies and believed to be partic- 
ularly suited near the Dirac point. 

The density of states within a narrow energy range 
near E/V — is shown in Fig.[7]for square samples of size 
L/a= 128, magnetic flux density B = 1/32 (h/e)a~ 2 and 
several disorder strengths /. Contrary to the zero mag- 
netic field case, the tails flatten with increasing disorder 
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FIG. 7: (Color online) The depression of the density of states 
within the central Landau band near energy E/V = 0.0 
for square samples of size L x = L y — 128 a and constant 
magnetic flux density p/q = l/32(h/e)a~ 2 . To enhance 
the differences, the DOS is shown on a logarithmic scale. 
The disorder strength for the various curves is f/(h/e) — 
0.01,0.02,0.03,0.05, 0.07, and 0.1, respectively. The topmost 
curve belongs to f /(h/e) = 0.01. 



for finite B, and the energy range of the DOS-depression 
gets broader. Please note that the DOS is plotted on a 
logarithmic scale for sake of clarity. A similar function as 
(H]) can be used to fit the energy dependence of the nar- 
row DOS depression, but now with a fitting parameter 
gb{f,L,B) which also depends on the constant magnetic 
flux density. A power-law relation gb(f,L,B) cx f~ 2 is 
obtained from the data shown in Fig. [7] This behav- 
ior is plotted in Fig. [H] for disorder values in the range 
0.01 < f/(h/e) < 0.1. The same dependence has been 
found also for size L/a = 192, magnetic flux p/q = 1/96 
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FIG. 8: (Color online) The disorder dependence of the fit- 
ting parameter gt(f,L,p/q). The size of the square sam- 
ples is L/a = 128 and the constant magnetic flux density 
B = 1/32 (h/e)a~ 2 . From the log-log plot, a power-law rela- 
tion gi,(f,L,p/q) = 355 (/ /(h/e))~ 2 is obtained. 
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FIG. 9: (Color online) The power-law dependence of the fit- 
ting parameter gi,(f,L,p/q) versus magnetic flux quanta p/q 
for square samples of size L/a = 192 (•) with f /(h/e) — 0.02 
and for L/a = 256 with disorder strength / /(h/e) = 0.03. 
The straight lines follow a cx (p/q) 1 ' 25 relation. 



and disorder strengths f/(h/e) =0.01, 0.02, and 0.05. 
Therefore, the energy range of the DOS depression broad- 
ens with increasing disorder strength and the tails be- 
come flat, which is completely opposite to the B = 
case. 

The specific magnetic field dependence of <?&(/, L, B) is 
not so easy to extract, because both the height and the 
width of the Landau band vary. In any case, the DOS 
tails become steeper with increasing B. For square sam- 
ples of size L/a = 256 and disorder strength / /(h/e) = 
0.03, a power-law relation gt(f,L,p/q) oc (1/g) 125 is ex- 
tracted from a fit to the relation corresponding to (g]) 
with q = 32,64, 128, and 256. A similar behavior is ob- 
served for L/a = 192, f/(h/e) = 0.02, and flux den- 
sities with q =48, 96, and 192. Both data fits can be 
seen in Fig. |]on a double logarithmic scale. This mag- 
netic field dependence means that the energy range of 
the DOS-depression becomes broader and therefore more 
important when the magnetic flux density gets smaller, 
approaching those applied in experiments. However, this 
power-law relation will probably only hold as long as the 
magnetic length Is = y/H/JeB) = a^J q/ (2np) remains 
smaller than the system size L. 

The size dependence of gt(f = 0.01 h/e, L,p/q = 1/32) 
is found to be cx L 2 5 in the range 64 < L/a < 192. 
Although the tails of the DOS-depression get steeper 
with increasing system size, p(E/V = 0) stays zero at 
the Dirac point in the range 64 < L < 512 investi- 
gated. Putting everything together, in the presence of 
a perpendicular magnetic field the dependence of the 
empirical fitting function on disorder strength, system 
size, and magnetic flux density can be summarized by 
g b Kf- 2 (L/l B f/ 2 . 
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FIG. 10: (Color online) The disappearance of the depres- 
sion in the density of states near the Dirac point due to an 
additional diagonal disorder potential of strength W/V — 
1CT 5 , 10~ 4 , 10~ 3 , 2x 10~ :! , and KT 2 , respectively. The system 
size is L/a = 128, the random flux disorder f/(h/e) — 0.05, 
and the constant magnetic flux density B = l/32(h/e)a~ 2 . 



IV. DISCUSSION AND CONCLUSIONS 

The density of states of non-interacting electrons mov- 
ing on a two-dimensional bricklayer lattice in the pres- 
ence of chiral symmetry preserving random flux disorder 
and a perpendicular magnetic field exhibits a narrow de- 
pression near the Dirac point at E/V = 0. The corre- 
sponding numerical results reveal a dependence on the 
disorder strength, the magnetic flux density, and on the 
size of the system. The latter is not simply a finite size 
effect because the special size dependence develops only 
near the Dirac point where p{E/V = 0) stays zero even 
though the steepness of the tails grows with increasing 
L. 

Since the DOS-depression can be removed by an ad- 
ditional diagonal disorder giving rise to inter-valley scat- 
tering, the origin of this feature must derive from the 
sub-lattice structure and the associated chiral symmetry 
of graphene's honeycomb lattice. The modeling of di- 
agonal disorder required an extra term J2 X y e x.yC x y c x,y 
in the Hamiltonian (Q]). The set of uncorrelated ran- 
dom disorder potentials {t x ,y\ was chosen to be box dis- 
tributed — W/2 < e XlV < W/2 with probability density 
1/W. The removal of the DOS-depression as a func- 
tion of additional short-range disorder potentials can be 
seen in Fig. [10] for systems of size L/a = 128, random 
flux disorder / /{h/e) = 0.05, and magnetic flux density 
B = l/32(h/e)a~ 2 . With increasing disorder strength 
W/V = 10- 5 ,10- 4 ,10- 3 ,2 x 10~ 3 , and 10~ 2 , the nar- 
row DOS-depression in the lowest (central) Landau band 
completely disappears. 

The occurrence of the DOS-depression in the absence 



of short-range diagonal disorder for B = could be ex- 
plained if the elastic scattering length diverged at the 
Dirac point. Then, for increasing system size a decreas- 
ing energy range around E/V = would exist where the 
elastic scattering length is larger than L. Within this 
energy interval, the transport would be almost ballistic 
and, due to the absence of scattering events, the DOS 
approaches the result of the ordered case and drops to 
zero at E/V — 0. If this size dependence of the DOS- 
depression were accessible in experiments, it would open 
the possibility for obtaining information about the elastic 
scattering length. 

Due to the lack of an analytical theory for the density 
of states of a disordered bricklayer model near the Dirac 
point, a relation proposed in the context of the thermal 
quantum Hall effec t 44 ' 45 and recently for disordered Dirac 
fermions on a honeycomb lattice^ was ventured. While 
the empirical relation ([4]) used to fit the numerical re- 
sults seems to work quite well, one has to keep in mind 
that only analytical calculations for a lattice model in 
the presence of random magnetic flux will eventually help 
to understand the complete situation. The usual way to 
start from the continuum Dirac equation may turn out to 
be not appropriate for a comprehensive description, if the 
depression in the density of states found in the present 
study would not show up in the former description. 

The implication of this observation and its impact on 
the scaling behavior and the critical properties at the 
Dirac point is evident, but still needs to be investigated. 
Usually, a noncritical density of states with a smooth 
energy dependence is assumed in the scaling analysis. 
In particular, the strong energy and size dependence re- 
ported above necessitate a reassessment of the conven- 
tional procedure applied in Ref. [l||. 

Although the disorder and magnetic field dependent 
depression found in the density of states at the Dirac 
point and the occurrence of the conductance peak split- 
ting reported previously 1 ^ are in agreement with several 
aspects observed in experiments mentioned in the in- 
troduction, it is clear that many-body effects and also 
single particle interactions like Zeeman splitting or spin- 
orbit scattering, which were not taken into account in the 
present investigations, may turn out to be the dominant 
effects in understanding these experiments. Nevertheless, 
the results of the calculations presented above may be 
helpful in finding out which type of disorder determines 
the electronic properties of real graphene samples. 
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